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1 Introduction 

Observations of precipitation over Greenland are limited. Direct precipitation measurements for 
the whole ice sheet are impractical, and those in the coastal region have substantial uncertainty 
(Bromwich and Robasky 1993) but may be correctable with some effort (Y ang et al. 1999). However, 
the analyzed wind, geopotential height and moisture fields are available for recent years, and the 
precipitation is retrievable from these fields by a dynamic method. Based on recent Greenland 
precipitation from dynamic studies (Bromwich et al. 1993; the National Centers for Environmental 
Prediction (NCEP)/ National Center for Atmospheric Research (NCAR) Reanalysis (Kalnay et al., 
1996); Chen et al., 1997a; European Centre for Medium-Range Weather Forecasts (ECMWF) 
Reanalysis for 1979-1993 (ERA- 15) (Gibson et al. 1997)), several deficiencies in the precipitation 
spatial distributions from these dynamic methods were evaluated by Bromwich et al. ( 1 998). Figures 
1 a and 1 b show the distribution of the mean annual Greenland precipitation for 1979-1 993 computed 
by the ERA- 15 and that for 1985-1999 retrieved recently by the co-equation method (Chen et al. 
1997a; Bromwich et al. 2001a), respectively. Figure lb is computed based on ECMWF operational 
data from TOGA Archive II of NCAR. 

The snow accumulation, M s , is determined by 

M = P - E - D -D (0 

iV1 s 1 sn mw 

where P and E are the rates of precipitation and evaporation/sublimation, and D sn and D mw are the 
horizontal mass divergence due to blowing snow and melting water, respectively. The precipitation 
is the most important source of accumulation, and on the other hand, the snow accumulation can 
be used to deduce the precipitation over ice sheets because it is the only quantity in (1) that can be 
measured at the ice core sites at the present time. Figure lc shows the mean annual precipitation 
distribution over Greenland based on the modern synthesis by B. Csatho et al. (2002, personal 
communication), who used corrected precipitation totals from the coastal weather stations in ad- 
dition to 253 accumulation values on the Greenland ice sheet. Box and Steffen (2001) showed that 
evaporation/sublimation is significant towards the coastal slopes. Their estimated evaporation/ 
sublimation values were added to the ice-sheet accumulation by B. Csatho et al. to approximate 
precipitation more closely in Fig. lc. Although there is a temporal discrepancy of precipitation 
between the long-term average in Fig. lc and that of recent years in Figs, la and lb, comparison 
of the modeled precipitation (in Figs, la and lb) with the directly and indirectly observed one (in 
Fig. lc) can still be used to check the accuracy of the modeled precipitation. 

It is easily seen that the distribution of precipitation in both Figs, la and lb is much smoother 
than that in Fig. lc. If the latitude 68" N is used to separate Greenland into its northern and southern 
portions, there are two relative large precipitation areas shown in Fig. lc in the western part of the 
northern Greenland. One is centered near the point (70" N, 47" W ) with a contour value of 50 an yr~ l 
and the other is centered near the point (76" N , 63" IV) along the north shore of Baffin Bay. Although 
a relatively large precipitation area located at the latitude 70" N can be found in both Figs, la and 
1 b, their central locations are a couple of degrees of longitude west from the relatively large center 
(at the point (70" A, 47" W)) as shown in Fig. lc. In both Figs, la and lb, there is no correspondingly 
large precipitation center near the point (76" N , 63" IV) along the north shore of Baffin Bay. In the 
east coast of the northern Greenland, Fig. lc has many mesoscale features, but both Figs, la and 
lb lack them. For example, there are several large precipitation centers located inland in the coastal 
region between 68" N and 78" /V and a slight one is located 16" N as shown in Fig. lc, but these 




relatively large areas cannot be found in both Figs, la and lb. In addition, there are also considerable 
differences in the western part of the southern Greenland between Fig. lc and both Figs, la and lb. 
The above comparisons show that the mesoscale features of the annual mean precipitation over 
Greenland in Figs, la and lb need to be further improved. 

The measured accumulation values from ice cores over Greenland have many unique features. 
The time scale of the measured accumulation is comparatively long. For example, only annual 
values can usually be obtained, thus the time scale of significant temporal variability is at least 
several years. For the annual time scale, the day-to-day variations are unimportant. However, the 
horizontal scale of the measured accumulation is very small. For example, the distance between 
some ice cores is a few kilometers, but the variations of the measured accumulation between them 
are often substantial. Thus, the horizontal scale of accumulation variations is only kilometers. 

The mesoscale characteristics of the Greenland topography is the major factor causing the 
mesoscale features of the precipitation shown in Fig. lc. During the computation of Fig. lb, a high 
(50 km grid length) resolution Greenland topography has been used (Bromwich et al. 2001a), but 
the retrieved precipitation still fails to capture the mesoscale features. This is because the data sets 
(TOGA Archive from NCAR or ERA data) of the analyzed wind, geopotential height and moisture 
fields used in Figs, la and lb must reflect the resolution of the topography adopted at the forecast 
center (here the ECMWF) from which the analyses originate; further the fields are sub-sampled to 
a resolution of 2.5" x 2.5". These large scale resolution (2.5" x 2.5") data sets are inconsistent with 
the high resolution Greenland topography used in the computation of Fig. lb. This inconsistency 
must have an important impact on producing the mesoscale features of the precipitation caused by 
the mesoscale topography. 

Anthes (1990) showed that a mesoscale model with realistic treatment of mesoscale topography, 
earth surface conditions and physical processes is capable of developing mesoscale precipitation 
and phenomena from large-scale initial conditions. Physically, this means that a mesoscale model 
can produce mesoscale systems after a certain integration time through interactions and feedbacks 
between the large-scale initial conditions and the mesoscale topography, earth surface conditions 
and physical processes. In the interaction processes, the large scale initial wind and temperature 
fields evolve to be consistent with the mesoscale topography and physical processes and become 
the ones having mesoscale features. 

The generalized co-equation method studied by Chen et al. (1997a) is very useful for retrieving 
precipitation to study the mass balance over ice sheets. A reasonable precipitation description on 
the annual time scale with a relatively high horizontal resolution can be obtained by this method 
much more easily than by global and limited-area models. However, the generalized co-equation is 
only a diagnostic relation; it derives precipitation immediately and does not include interactions 
and feedbacks between the large scale initial conditions and mesoscale topography even if the 
mesoscale topography is correctly specified. As shown in Fig. lb, this method does not have the 
ability to generate enough mesoscale features of the Greenland precipitation from the large-scale 
initial conditions, which is an weakness of the co-equation method in comparison to mesoscale 
models. In this paper, this diagnostic method is further developed and improved by using an iterative 
method to adjust the large-scale 2.5" x 2.5" resolution analyzed data to be partially consistent with 
the high resolution Greenland topography. A balanced divergence equation is solved with the fixed 
external wind boundary condition (Chen et al. 1996) by iteration, and the new method is referred 



to as the iterative balanced divergence equation method. 

In a mesoscale model, the precipitation and atmospheric motion over mountainous regions, 
especially near steep slopes of mountains and ice sheets, are greatly influenced by the computational 
accuracy of the horizontal pressure gradient force over these regions. Colie et al. (1999) evaluated 
the 36- and 12-km resolution Penn State/NCAR mesoscale model (MM5) (Grell et al. 1994) pre- 
cipitation forecasts and NCEP’s 10-km resolution Eta Model (Eta- 10) forecasts across the Pacific 
Northwest of the U.S. and found that the 12-km MM5 tends to generate too much precipitation 
along the steep windward slopes. The Eta- 10 overpredicts precipitation along the windward slopes 
more than the 12-km MM5 even though the step-mountain approach of the Eta-coordinate system 
is used. The physical parameterizations of the MM5 have been adapted for applications over polar 
ice sheets (Bromwich et al. 2001b; Cassano et al. 2001) and the modified code is termed Polar 
MM5. Cassano et al. (2001) used the Polar MM5 to simulate a complete annual cycle from April 
1 997 through March 1 998 over the Greenland Ice Sheet. The simulations are a series of 48-h forecasts 
of which only the forecasts of the second 24-h period are used to represent each day. The modeled 
precipitation distribution is excessive along the steep coastal margins with spot values in excess of 
400 cm yr~ x located on the southeast coast while the corresponding observed values are close to 
120 cm yr~ l . The precipitation errors are similar to those over the Pacific Northwest found by Colle 
etal. (1999). Chen and Bromwich (1999) (hereafter referred to as CB99) showed that the horizontal 
pressure gradient force in sigma-coordinates can be computed more accurately by separating this 
horizontal vector into its irrotational and rotational parts, which are expressed by the equivalent 
geopotential and geo-streamfunction, respectively. This method has been successfully used by Chen 
et al. ( 1 997a) and Bromwich et al. ( 1 999; 200 1 a) in their (o-equation method to compute precipitation 
over the steep slopes of the Greenland ice sheet. Recently, this method has also been utilized in the 
MM5 to improve precipitation prediction over the steep slopes of mountains and ice sheets (Chen 
et al. 2001a, b). 

In Section 2, the iterative balanced divergence equation method to adjust the large-scale 2.5" x 2.5" 
resolution analyzed data to be partially consistent with the high resolution Greenland topography 
is discussed, and a test of an annual precipitation over Greenland is shown by the method with and 
without iterations. In order to test the capability of this method for retrieval of the mesoscale features 
of the precipitation affected by topography over polar ice sheets, the annual mean precipitation 
distributions over Greenland from 1985-1999 computed by this method are shown in this section. 

Because rises and falls of the ice surface elevation over a relatively short period are primarily 
produced by the snow accumulation, comparison of the interannual trend of precipitation with that 
of the observed surface elevation change over the Greenland ice sheet is discussed in Section 3. 

2 An iterative balanced divergence equation method for improved retrieval of the meso- 
scale features of the annual precipitation over Greenland 

a. The equivalent geopotential and geo-streamfunction in a coordinates 

The horizontal wind can be separated into its irrotational and rotational parts and expressed by 
V = -V% - k x Vy 



where % and \|/ denote velocity potential and streamfunction, respectively. 

The vertical coordinate a is defined by G = pip*, where p,(x,y,t) is the surface pressure. The 
horizontal pressure gradient force (HPGF) G in G-coordinates is expressed by 

G = -V<\>(x,y,o,t)-RTVlnp*(x,y,t) (2) 

where §(x,y,o,t) is the geopotential in G coordinates. Here G is also a horizontal vector. CB99 
proposed that the vector G can also be separated into its irrotational and rotational components and 
expressed by 

G = -V^ - k x Vr| (3) 

where and r|(jc,y,G,0 are referred to as an equivalent geopotential and a geo- 

streamfunction, respectively. 

Because the divergence of the HPGF, V ■ G, from the right hand side of (3) must be equal to that 
from the right hand side of (2), we have 

. d L„, a/np.yaf, _ dlnpA 


V‘ty e (x,y,C,t) = V‘ty(x,y,0,t) + ^- RT(x,y,o,t ) 


+ 5 - RT(x,y,a,t)- 


In order to reduce computational errors in a limited region, a harmonic-sine spectral method 
(Chen and Kuo 1992) is used. Thus, the geopotential and equivalent geopotential in G-coordinates 
can be separated into their inner and harmonic parts as 

<|> = <t >,(*,>’, o, t) + <j ) h (x,y, o, t), <|> e = <t > ei (x,y,o, t ) + <t) f/l (x,y ,G, t) (5) 

where the inner part of the equivalent isobaric geopotential height in G-coordinates, (j)^., can be 

derived based on (4) from the solution of the following Poisson equation 


V> ei = V3|),+— RT 


with zero Dirichlet boundary value. 

Both -V^Cx.y.G.r) and -V$(x,y,p,t) are the irrotational part of the HPGF but they are in o- 
and p-coordinates, respectively. Thus, the equivalent geopotential ^(x.y.o, t) can be used in 
G-coordinates in the same way as <ty(x,y,p,t) is used in p-coordinates. The divergent and rotational 
parts of the HPGF are expressed by V 2 ^ and V 2 r|, respectively, and the divergence of the Coriolis 
force is denoted by f Q £l, where / = / 0 +/, where f 0 is the average value of the Coriolis parameter in 
the integration region, and/* is its deviation from^- Here Q. is vorticity denoted as (2.1) of CB99. 
The divergence of the HPGF can be separated into its geostrophic and ageostrophic parts, where 
the geostrophic part is equal to f 0 Q as 

V\ g =f„PL 

and the ageostrophic part of V 2 ()) t , is the difference between V 2 ^, and its geostrophic part and denoted 
by 

v\, ( = v 2 4> t . -/ 0 n 

Thus the divergence of the HPGF is the sum of its geostrophic and ageostrophic parts as 
= V\ g + 


(7) 



In the quasi-geostrophic approximation, the ageostrophic part of the divergence of the HPGF, V 2 <j) f 
always vanishes, thus we have 

v 2 i = v^, # =/^=/„vV {8) 

b. The vertical finite difference forms of the vorticity, divergence, continuity, hydrostatic and 
thermodynamic equations in 0 coordinates 

The vertical distribution of the variables is shown in Fig. 6 of CB99. In the present paper, 22 
0-levels at 0=0.010, 0.030, 0.050, 0.070, 0.090, 0.1 10, 0.140, 0.180, 0.220, 0.270, 0.340, 0.420, 
0.505, 0.590, 0.670, 0.750, 0.820, 0.875, 0.920, 0.955, 0.980, and 0.995 are used in the vertical. 
Based on (5.9) and (5. 10) of CB99, the vorticity and divergence equations are expressed by 

dQ, X — i i (9) 


a t 

dD i 

dt 


= -foD X +£l adr 4 
=f 0 nX-v\X-v 2 EX+D adv X 


( 10 ) 


where the column vector X X of a variable X is denoted by 

Xl=(X„ .... X k , ..., X N f (11) 

and ( . . . ) T is for transpose, and N is the total number of 0 -levels. Here the term E, X is denoted by 
£. X= [i m\U 2 + V 2 ),/2] X, and the terms Q a(lv I and D aih , X are the variation rates of the vorticity and 
divergence caused by advection and denoted by (5.6) and (5.7) of CB99, respectively. 

Based on (6.1) of CB99, the vertical finite difference form of the continuity equation is written 
as 

3 In p * 


dt 


+ mlUD 1= P 


adv 


where D X is the column vector of horizontal divergence and IT indicates the row vector as 


n = (A 0 „ 


Ao, 


A0„). 


( 12 ) 


(13) 


The map scale factor is separated by m~ = (m‘) 0 + (nr)', and (tn~) 0 and (m~) are the average value 
in the integration region and deviation, respectively. The term P adv of (12) is dominated by the 
surface pressure advection and it is expressed by 


P = 

r adv 
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(14) 


dx J dy 

The finite difference form of the hydrostatic equation based on (6.7) of CB99 is expressed by 

4» != <>. X +RBT X (15) 

where matrix B is a upper-triangular matrix shown by (A.29) of Chen et al. (1997b), and <]>» X- <j).I. 
Here I is the identity matrix and (j). = gH », where //, is the height of the earth’s surface. 

The temperature can be separated into two parts 
T( x ,y,a,t) = T o (0) + T\x ,y,o,t) = T 0 (o) + ( T(x , y , 0, /) - T 0 (a)) 


( 16 ) 


where T 0 {<3) is the averaged value of T(x, y, o, t) over the constant <5 — surface, and 7 ' is its deviation. 
The temperature separation (16) is different from that (5.2) and (6.9) of CB99 but is the same as 
that used by Chen et al. (1997b). The separation (16) is more natural than separation (5.2) and (6.9) 
used by CB99. Based on (2.28) of Chen et al. (1997b), the vertical difference form of the ther- 


modynamic equation is written as 


where 


T ha A=T ad A-(m 2 y¥Di+P T i 


(17) 

(18) 


Here the matrix F is denoted by (A. 34) of Chen et al. (1997b), and P T is the diabatic heating. 


c. The equations of the equivalent geopotential and ageostrophic part of the HPGF divergence 

Substituting the inner part of the hydrostatic equation (15) into (6), the vertical difference form 
of (6) is expressed by 


V 2 ^ 1= v 2 <), i +RBV 2 T i +^-( RT ' 

ox 1 dx 
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Taking the partial derivative of (19) with respect to t, and utilizing (16), we have 
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( 20 ) 


Comparing the above equation (20) with the similar equation (6.16) of CB99, there are additional 
two more last terms on the right hand side of (20). Substituting (12), (15), (16) and (17) into (20), 
then the equation of the HPGF divergence (20) becomes 

av 2 4> t . I 


dt 


■ + m~AD !- V-<P eMd i 


( 21 ) 


where matrix A is 


A = R [B¥ + T 0 i n) (22) 

and V 2 <J> t , ha(l i denotes the variation rate of the HPGF divergence caused by the advection and 
heating. For a variable X in a limited region, we have 

X = X,+X h , (23), V 2 X, = 0. (24) 

Thus, V 2 d> t , w i can be rewritten as 

4.= V 3 <t».., w , i= /?BV 2 7 W , l+RT 0 l V 2 />„,, 
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(25) 

dx dx dy dy 

There are two differences between the above equation (25) and similar equation (6.19) of CB99. 
One is that the variation rate of the HPGF divergence is used in (25) rather than that of the inner 



part of the equivalent geopotential. The second is that there are two more additional terms on the 
right hand side of (25) than on that of (6. 1 9) in CB99. With these two additional terms, the variation 
rate of the HPGF divergence caused by the advection and heating, V'<t > c w 1, can be described 
more accurately. 

Based on (7), the equation of the ageostrophic part of the HPGF divergence can be derived 
from (9) and (2 1 ) and expressed by 


av 2 <j) f , fl i 

dt 


+ AV'D -l -f 0 D si V ^ haj a ■i 


(26) 


where 


V 2 <£> . , 1= V 2 <t> . . I 

e,haci,a e,had 




(27) 


is referred to as the variation rate of the ageostrophic part of the HPGF divergence caused by 
advection and diabatic heating. 


d. The balanced divergence equation in (5-coordinates 

If the tendencies of the divergence and ageostrophic part of the HPGF divergence in (10) and 
(26) are neglected, this approximation is referred to as a balanced ageostrophic approximation (Chen 
et al. 1996). In this case, equation (26) becomes 

m 2 \V 2 D t si -f(Pi sl= V 2 <tv w>a i ■ ( 28 ) 

Equation (28) is referred to as the balanced divergence equation in o -coordinates, and it is also a 
divergence form of the generalized co-equation. In this equation, the diabatic and advection terms 
computed by the ageostrophic wind are the same as those in the generalized co— equation in p- 
coordinates (Pauley and Nieman 1992), but the effect of orography on the vertical motion is much 
better described than that in p-coordinates. 

Equation (28) can be transformed into the equation of its vertical mode. We introduce a matrix 
E in order that the following relation is satisfied 

E~' AE = G = diag(gh gh 2 , ... ,gh N ) (29) 

Let the velocity potential in the vertical mode and physical space are expressed by 
D,, -1= E“'D, I, and D, 1= ED,, l (30) 

respectively. Equation (28) is multiplied from left by the matrix E“', and then its equation of the 
vertical mode separately is written as 


Cy~D t , k -fiDn = V 2 d> t , , ia , +f~D htk , 

Xr 

II 

NJ 

, N 

(31) 

where 







(32) 

and Q. = m 0 sjgh k 



(33) 


Here C k is the gravity-inertia wave phase speed for the k-th vertical mode, and 

2 ,n o8 h k ( £* Y 

"" ~UJ 


(34) 



where L^ k is the radius of deformation of the k-th vertical mode. 

From the solution D, of (28), the divergence is computed by D = D, + D h . Using the continuity 
equation and vertical finite differencing, the pressure vertical velocity © in a coordinates is expressed 
by 


. , ( . 3lnp» . 3lnp» 

l=m 2 (I-C) U i^r~+V 

l ox oy 


— m 2 CD X 


(35) 


where C is a lower-triangular matrix and shown by (A. 15) of Chen et al. (1997b). 

The condensation method and the procedure for computing precipitation from the vertical motion 
are the same as those presented by Bromwich et al. (2001a). 


e. An iterative solution of the balanced divergence equation with the fixed external wind boundary 
value 


As pointed out in section 1 , the ©-equation method cannot generate enough mesoscale features 
of the Greenland precipitation from the large-scale (2.5° x 2.5") data sets even if the mesoscale 
topography is specified. This is because the large scale resolution data sets are generated from the 
same large scale topography and are not consistent with the high resolution topography. During the 
interpolation from the large scale (2.5" x 2.5") data in p-coordinates to the high resolution ones in 
a-coordinates, the surface pressure and temperature are computed based on the given high resolution 
topography and hydrostatic equation, and they are more consistent with the high resolution to- 
pography, especially the surface pressure. The horizontal wind is directly interpolated from the 
large scale p-coordinates to the new high resolution a-coordinates, and it may not be consistent 
with the mass field (temperature and surface pressure) in a-coordinates based on the high resolution 
topography. The balanced divergence equation (28) or (31) is one of the balanced equations used 
in the implicit nonlinear normal mode initialization (NNMI) (Temperton 1988) or balanced 
ageostrophic initialization (Chen et al. 1966) except that the equivalent geopotential and a-coor- 
dinates are used here. The NNMI or balanced ageostrophic initialization is through an iterative 
method to adjust the unbalanced ageostrophic wind to be the balanced ageostrophic wind, and the 
initialized wind is the sum of the geostrophic and balanced ageostrophic winds (Chen et al. 1996). 
This basic result will not change if the equivalent geopotential and a-coordinates are used. In the 
©-equation method (Chen et al. 1997a; CB99), the vertical motion © is computed from (35), in 
which only the divergence D is adjusted once through the solution of (3 1 ) but the wind components 
used in the advection term do not change. Thus, the vertical motion ©might be not adjusted enough 
and not consistent with the high resolution topography in this method. In the present section, the 
solution of the balanced divergence equation (31) is further adjusted by an iterative method as 


follows. 

For the vth iteration, equation (3 1 ) can be written in the form 
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(36) 


The harmonic part of the divergence Dff in (36) is assumed to be unchanged in the iteration. The 
nonlinear term, is calculated by using the values of U iv ~ 0 X, F (v 0 i and D <v “ " i at 

the (v - l)th step. At the beginning step, v = 0. 


Equation (36) is solved easily by the double sine series due to their homogeneous boundary 
value for the inner parts, and its solution is directly expressed by 
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where F _l [ - ] is the inverse double sine Fourier transform operator denoted by (A. 3) of CB99. 
During the above iteration, only the divergent component of the wind is modified within the region. 
Because the external wind is both nondivergent and irrotational in a limited region, no matter how 
the divergence is modified within the region, the external wind does not vary in the region and up 
to the boundary (Kuo and Chen 1992; Chen et al. 1996). Thus, it is very natural to use a fixed 
external wind as the lateral boundary condition in the iteration. If the interpolated wind components 
in a limited region bounded by L are given by the column vector as 

U i0 \x,y)i, V i0 \x,y)i. (38) 

The external wind can be derived by the natural method (Chen et al. 1996) as follows. The vorticity, 
£2 <0) 4, and divergence, D m 4, are computed from the wind components (39). The inner parts of the 
streamfunction and velocity potential are obtained by solving the Poisson equations 

V>; 0) 4=f2 (0) 4, V 2 x‘ 0) 4= D (0) 4 (39) 

with the homogeneous Dirichlet boundary values. The internal wind is computed by 


a\l/ (0) 4 r>/ 0) 4 
up 4=- : .- + L 


avi/ <0, 4 3 y (0) 4 
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(40) 


dy dx 
and the external wind is derived by 

uf> 4= u m 4 -up 4, v? 4= v (0) 4 -vf ] 4 . (4 1 ) 

In each iteration, it is necessary to compute the nonlinear advection terms, Q. ady and Mv in 
(27). These nonlinear terms are computed by the transform method (Orszag 1970; Eliasen et al. 
1970). 

At v-th step, the inner part of the divergence £>/ v) is derived. The divergence is derived by 


D {v) = D- v) + D 


( 0 ) 


Thus, the inner part of the velocity potential can be derived from the Poisson equation 
v 2 x< v) = d , v) . 

with zero boundary value. The internal wind at v-th step is computed by 

L «n 
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The total wind is reconstructed by 

u 4 <v) = uf 4 +up 4, v /(v) 4= vf 4 +v?' 4 

where the external wind components. Up and V^ 0 ’, are obtained from (41). 
In the iteration, equation (35) can be expressed by 
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f An example of the annual precipitation for 1986 over Greenland 


(42) 


(43) 
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(45) 


(46) 
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In order to check how the mesoscale features of the computed precipitation is affected by the 
iteration, the annual precipitation for 1986 over Greenland is computed with and without iteration, 
respectively, but the other parameters are all the same. The computational domain and topography 
of Greenland are presented in Fig. 2a, and the modern terrain dataset of Ekholm (1996) for Greenland 
topography is used. The mesh size is 1 1 1 x71, and grid spacing is 25 km. 

There are three methods which can be used to compute the pressure vertical velocity to based 
on (35). The first method is to compute the pressure vertical velocity 0) directly by 

'^)l= m\ I - C) { U {0) i + V (0) i ^ 1 - m 2 CD m i (47) 

\ p ) \ ° x °y J 

where (7 <0) i and F ,0) i are interpolated from the 2.5" x 2.5" resolution data in p-coordinates to the 
25 km resolution grid of the limited region in o-coordinates, and the divergence D m i is computed 
from t/ (0) i and T (0) i. The method (47) is referred to as a kinematic method. The horizontal di- 
vergence derived from this method is very sensitive to small errors of the interpolated wind com- 
ponents f/ (0) 4- and F (0) 4-. 

The second method is to derive pressure vertical velocity by 
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where the divergence D (i) 4- is derived from the divergence equation (28), in which the forcing term 
V 2 <U e had at 4- based on (27) has two terms of similar magnitude. Thus, the divergence D 0) 4- derived 
from (28) is much less sensitive to the wind errors than D (0) 4- computed from the kinematic method. 
The method (48) may be referred to as the co-equation method without iteration, and it was used by 
Chen (1997a) and Bromwich et al. (1999, 2001a) except that D (,) i is derived from the velocity 
potential form of the equation (28). In this method, the vertical velocity co (l) i is more accurate and 
more consistent with the high resolution topography than co <0> 4-. 

In the third method, the divergence D u> i is further used to improve divergent component of 
the wind by using (44) and (45) to obtain U 0) 1 and V 0) 1. Then f/ (1) i and V (1) i can be used to 
improve the horizontal advection in the forcing term of (35) to derive D (2) 4-, and the vertical velocity 
cd < 2) i is computed by 
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(49) 


The third method based on (49 ) is referred to as an iterative balanced divergence equation method. 
In general, this iterative method converges rapidly, and it is accurate enough to let the value of v 
be 2 or 3. In the present paper, that v=2 is used in the iterative method. 

As an example, the distribution of the annual precipitation for 1986 over Greenland computed 
from (48) with the 25-km grid length without the iteration is shown in Fig. 2b. Comparing Fig. 2b 
with Fig. lb, it is seen that, although the grid length used in Fig. 2b is a half of that used in Fig. lb, 
the distribution of precipitation in Fig. 2b is quite similar to that in Fig. lb and its mesoscale features 
are not improved. Figure 2c depicts the annual precipitation of the same year computed by the 
iterative balanced divergence equation method based on (49). Figure 2c has many mesoscale features 
of the precipitation, especially a relatively large area near the point (70" N, 47" VT) and a center near 
the point (76" N, 63" W) along the north shore of Baffin Bay in the west coastal region of Greenland. 
In the east coast of the northern Greenland, there are two relatively large mesoscale precipitation 
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centers located inland in the coastal region between 68" TV and 78" TV’ and similar to Fig. lc, the 
smaller one is located 76" TV as shown in Fig. 2c. Thus, the precipitation in Fig. lc is more similar 
to the annual precipitation in Fig. 2c than to that in Fig. 2b. This example shows that the iterative 
balanced divergence equation method has capability of improving the mesoscale features of pre- 
cipitation over Greenland affected by the high resolution topography. 

g. The retrieved mean annual precipitation distribution of the Greenland for 1985-1999 

The mean annual precipitation distribution over Greenland computed by this method for 1 985-99 
is shown in Fig. 3, which is calculated based on the ECMWF operational data (TOGA Archive II 
from NCAR), the topography of Fig. 2a at 25 km resolution and the 22 levels in the vertical of the 
model. 

Comparing Fig. 3 with Fig. lb, it is seen that the distribution of the mean annual precipitation 
in Fig. 3 has much more mesoscale features than that in Fig. lb. The two relatively large accumulation 
areas, which are centered near the point (70" TV, 47" VF) and along the north shore of Baffin Bay at 
the point (76" TV, 63" IV), respectively, in the western part of the northern Greenland shown in Fig. 
lc, have been retrieved in Fig. 3, but they are not simulated well by the method without iteration 
as discussed in Section 1 . The areas encircled by the contours of 10 and 20 cm y~' over the central 
Summit region depicted in Fig. lc are more similar to those in Fig. 3 than in Fig. lb. For example, 
the area around by the closed contour of 10 cm y _l in Fig. lb is too large, which means that the 
modeled result in Fig. Ibis too dry over the central Summit region. This weakness has been improved 
in Fig. 3. In the east coast of the northern Greenland, there are two relatively large mesoscale 
precipitation centers located inland in the coastal region between 68" TV and 78" TV in Fig. 3: the 
strong one is located at 76" TV while the weak one is located at 76" TV. These two mesoscale areas 
are similar to those in Fig. lc, but they are not retrieved in Fig. lb. In addition, the significant 
precipitation errors in the south of 65" TV of the west coastal region of the southern Greenland in 
Fig. lb (in comparison with Fig. lc) have also been corrected in Fig. 3. The mesoscale distributions 
of both Figs. 3 and lc in this region are quite similar. However, the differences of the mesoscale 
distributions along the west coast between 65" TV and 72" TV in Figs, lc and 3 are different, and these 
mesoscale differences of the precipitation in this region need to be further studied. From all of the 
above comparisons, it is seen that many features of the annual mean precipitation over Greenland 
which could not be retrieved by the method without iteration as shown in Fig. lb have been greatly 
improved in Fig. 3. 

3 Comparison of the interannual trend of the precipitation with the observed surface 
elevation change for 1993-1999 over the Greenland ice sheet 

a. The equation of variation of the ice sheet surface elevation 

The variation rate of the surface elevation of an ice sheet, dH i( Jdt, is determined by 

C ^ = M-D f (50) 

dt ' 

where M s , and D f denote the variation rates of the surface elevation caused by accumulated snow 
and firn densification, respectively. Here M s is expressed in snow depth rather than in water 
equivalent. The firn densification only changes the snow surface elevation but does not alter the 



net water equivalence, and it may also be referred to as snow deflation. Based on ( 1), equation (50) 
can be rewritten for a relatively short time period as 


dH ue 

dt 


= P - E - D m - D -D 


(51) 


As D s in (50), the terms P, E, D„ and D mw in (51) are all expressed in snow depth rather than in 
water equivalent. Because the precipitation is a very important source in accumulation, it is also 
very important for the change of the surface elevation of ice sheets. In order to understand what is 
responsible for the changes of the surface elevation of the Greenland ice sheet, it is necessary to 
study the corresponding precipitation changes. 

b. The interannual trend for 1993-1999 of the Greenland precipitation at elevations above 2000 


meters 

Recent advances in airborne laser altimetry and global positioning system (GPS) technology 
have made possible the large-scale assessment of elevation change characteristics of the entire ice 
sheet through repeated surveys separated in time. Such repeated surveys in 1993 and 1998 (Krabill 
et al. 1999) showed that the southeast margin of the Greenland ice sheet has been thinning. Aircraft 
laser-altimeter surveys over northern Greenland in 1 994 and 1999 have also been studied by Krabill 
et al. (2000), and they reported changes in the surface elevation of Greenland between 1993 and 

1999 derived from radar and laser altimetry and estimated coastal melting. It is found that, above 

2000 meters elevation, the entire ice sheet is in balance on average but has some regions of local 
thickening or thinning. The changes of surface elevation of the Greenland ice sheet above 2000 m 
are shown in Fig. 4a. 

Above 2000 m surface elevation, most of the northern ice sheet lies above the region of summer 
melting. McConnell et al. (2000) derived changes of the ice-sheet elevation in southern Greenland 
for the years 1978-88, using a physically based model of firn densification and records of annual 
snow accumulation reconstructed from 1 2 ice cores at elevations above 2000 meters elevation. They 
found that the patterns of elevation change derived from snow accumulation agree closely with 
contemporaneous satellite measurements of the surface elevation change of the ice sheet. Thus, the 
effects of melting on the surface elevation change of southern Greenland above 2000 m should also 
be small. 

In order to compare the temporal variability of precipitation with that of the surface elevation 
over Greenland above 2000 m, the spatial distribution of the slope of the linear regression line of 
the annual precipitation from the balanced divergence method for 1993-1999 has been computed. 
A color-coded figure of the slope of the linear regression line of the annual precipitation from 
1993-1999 over Greenland above 2000 m is shown in Fig. 4b. It should be pointed out that the unit 
used in precipitation is cm/year in water equivalent, and a multiplying factor of about 3.3 is necessary 
to transform the water equivalent values to the thickness of snow (R. Thomas, personal commu- 
nication, 2001) in Fig. 4b. 

In southern Greenland (south of 70" A0 above 2000 m, there are three thinning areas in both 
Figs. 4a and 4b. Two thinning regions are located to the south of 61" N over its eastern part and 
western part, respectively. The strength and area of these thinning regions over the eastern part are 
larger than those over the western part. It should also be noted that the colors of contour spacing 
used in Figs. 4a and 4b are not the same. There are two other thinning regions centered in 
69.5" N, 47" W and 69.5" ^,34" W, respectively shown in both Figs. 4a and 4b. There are three 
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Figure 4 a) A part of Fig. 4c for Greenland ice-surface elevation change in the region above 
2000 m; b) A part of Fig. 4d for annual precipitation trend in the region above 2000 m; c) 
Greenland ice-surface elevation change dh/dt in unit cm/year for 1993-1999 derived mostly 
from airborne laser-altimetry. The 13 coastal stations are shown in green along with the dh/dt 
value derived from the PDD anomalies (Krabill et al. 2000); d) the annual precipitation trend 
over Greenland computed by the iterative balanced divergence method in unit cm/year. 






thickening areas located at about (68" N,3T' (65"/V,47" Vy) antl (62.5" N, 48" VP) s h° wn i n 

Figs. 4a and 4b. Thus, the central locations of thinning and thickening areas in the southern Greenland 
above 2000m in Figs. 4a and 4b are in good agreement with each other. 

In northern Greenland (north of 70" (V) above 2000 m, there are two major thickening regions 
centered at about (72" N, 48" W) and (76" N, 28" W) and two weak thinning areas centered at about 
(76" N, 44" IV) and (73" N, 29" W) shown in Fig. 4a. These two thickening and two thinning regions 
correspond to two positive and two negative areas of precipitation change respectively centered at 
about the same locations in Fig. 4b. One positive precipitation increase region centered at about 
(71" AT, 32" W) shown in Fig. 4b is not matched to a thickening area in Fig. 4a. This may be due to 
deficiencies in the interannual variation of the precipitation retrieval or other causes, for example, 
relatively large melting in this region. 

From the above, it is seen that the Greenland ice sheet above 2000 meters elevation from 
1993-1999 is nearly in balance with some regions of local thickening or thinning, and the 
altimetry-derived ice-sheet thickening and thinning are approximately consistent with the precip- 
itation change. This situation was first discussed by Bromwich et al. (2001a). Because some me- 
soscale features of the precipitation over Greenland have been improved by the balanced divergence 
equation method, the newly computed annual precipitation trend is compared again with the 
observed changes of the ice sheet surface elevation, and the results are in better agreement with the 
measured local thinning and thickening. 

c. The interannual precipitation trend for 1993-1999 at the elevation below 2000 meters of 

Greenland and the downward trend for 1993-1999 over southern Greenland 

At elevations below 1700 m, radar altimeter data become unreliable (Thomas etal. 1999). Krabill 
et al. (2000) calculated a hypothetical thinning rate at the coast on the basis of the coast positive 
degree day (PDD) anomalies, using a factor of 9 mm per PDD. From their approach, only melt is 
considered near coast in the thinning rate. They excluded measured surface elevation change in 
coastal areas; instead, an interpolation was used between the calculated PDD thinning rates due to 
melting and nearest observed elevation changes to yield thinning rates over the ice-covered coastal 
regions. Based on their method, the annual trend of the surface elevation of the Greenland ice sheet 
including the region above 2000 m is shown in Fig. 4c. Below 2000 m surface elevation, Fig. 4c 
shows that thinning predominates at lower elevations along about 70% of the coastal regions, with 
rates about 1 meter per year close to the coast. The thinning rates exceeding 1 m/year over the 
coastal areas shown in Fig. 4c are probably too large to be caused by melting only. 

Based on (5 1 ), changes of surface elevation over the ice sheet are due to not only the melting 
but also precipitation, evaporation/sublimation, firn densification and drifting snow. Thus, it is not 
appropriate that only melting is considered in estimating the change of surface elevation of ice sheet 
even in the coastal areas. In order to get a relatively accurate estimate of surface elevation change 
over the whole Greenland ice sheet, especially the region below 2000 m, not only precipitation but 
also the melting, evaporation/sublimation, firn densification and drifting snow need to be studied 
in the modeling and diagnostic studies. 

The slope of the linear regression line of the annual precipitation from 1 993- 1 999 over Greenland 
including the regions above and below 2000 m is shown in Fig. 4d. It is clear that the linear trends 
over the coastal region on the average is negative, especially with relatively large negative values 
in the coastal regions of southern Greenland. The linear trends in annual precipitation from 



1985-1999 for Greenland have also been calculated, similar negative values are also found in the 
coastal regions of southern Greenland (figure omitted). Thus, there is a significant downward trend 
in annual precipitation from 1985-1999 for the southern Greenland and its coastal regions. This 
result is consistent with the report of Bromwich et al. (1999) that a significant downward trend in 
annual precipitation from 1985-1995 for all of Greenland and its southern and central-west coastal 
regions, amounting to 3% per year had been retrieved by the dynamic method (Chen et al. 1997a). 
The similar results for 1985-1 999 were also shown by Bromwich et al. (200 1 a) based on the improved 
dynamic method. 

4 Conclusion 

Based on the studies shown in the above sections, the following conclusions can be reached. 

(1) The advatage of the generalized co-equation method in o-coordinates developed by Chen et 
al. (1997a) and CB99 is to obtain a reasonable precipitation over steep slopes of the Greenland ice 
sheet in the annual time scale more easily than global and limited-area models, but its weakness is 
not able to generate the mesoscale features of the Greenland precipitation very well from the 
large-scale initial conditions. The generalized co-equation method is further developed now by using 
an iterative method to adjust the large-scale 2.5° x 2. 5" resolution analyzed wind to be partially 
consistent with the high resolution topography. In this method the balanced divergence equation is 
solved by iterations with the fixed external wind lateral boundary condition, thus, this method is 
further referred to as the iterative balanced divergence method. The computed results show that this 
iterative method has good capability in computing the mesoscale features of the annual mean 
precipitation affected by high resolution topography from the large scale analyzed wind. 

Many mesoscale features of the mean annual precipitation distribution over Greenland are 
retrieved by the iterative balanced divergence method. For example, the two relative large pre- 
cipitation areas, which are centered near the point (70" N, 47" W) and near the point (76" N , 63" VF) 
along the north shore of Baffin Bay in the western part of Greenland, are retrieved, but they are not 
simulated by the original co-equation method without iteration. The computed mean precipitation 
over the central Summit region is also improved to be quite similar to the observed precipitation. 

(2) The Greenland ice sheet above 2000 meters elevation from 1993-1999 is nearly in balance 
with some regions of local thickening or thinning, and the altimetry-derived ice-sheet thickening 
and thinning are approximately consistent with the precipitation change. Because some mesoscale 
features of the precipitation over Greenland have been improved by the balanced divergence 
equation method, the computed annual precipitation trend is in good agreement with the measured 
local thinning and thickening. In order to get a relatively accurate estimate of surface elevation 
change over the whole Greenland ice sheet, especially the region below 2000 m, not only precip- 
itation but also the melting, evaporation/sublimation, firn densification and drifting snow need to 
be studied in the modeling and diagnostic studies. 

The linear precipitation trend over the coastal region on the average is negative, especially with 
a relatively large negative value in the coastal regions of southern Greenland. There is a significant 
downward trend in annual precipitation from 1985-1999 for the southern Greenland and its coastal 
regions. 
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